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Abstract 

An algorithm based on the Ehrhch-Aberth root-finding method is presented for 
the computation of the eigenvalues of a T-palindromic matrix polynomial. A 
structured linearization of the polynomial represented in the Dickson basis is 
introduced in order to exploit the symmetry of the roots by halving the total 
number of the required approximations. The rank structure properties of the 
linearization allow the design of a fast and numerically robust implementation of 
the root-finding iteration. Numerical experiments that confirm the effectiveness 
and the robustness of the approach are provided. 
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1. Introduction 

The design of effective numerical methods for solving structured generalized 
eigenvalue problems has recently attracted a great deal of attention. Palin- 
dromic matrix polynomials arise in many applications [1^. An n x n matrix 
polynomial of degree k P{z) = Y!1^qA,z\ Ak ^ 0, A, e C"^", < i < fc^s 
said to be T-palindromic if Aj — Ak^i for i = 0, . . . , k. It is well-known 



21| that the palindromic structure induces certain spectral symmetries: in par- 
ticular if A ^ is an eigenvalue of P{z) then 1/A is also an eigenvalue of P{z). 
Numerical solution methods are generally asked to preserve these symmetries. 

The customary approach for polynomial eigenproblems consists in two steps: 
First P{z) is linearized into a matrix pencil L{z) = zX + Y, X,Y e £nkxnk ^ 
and then the eigenvalues of L{z) are computed by some iterative solver. The 
usual choice of the matrix QZ algorithm applied to a companion linearization 



12| of P{z) is implemented in the Matlab function polyeig. An alternative solver 
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based on the Ehrlich-Aberth root finding algorithm is proposed in [5j for dealing 
with certain structured linearizations. Specifically the method of [5| is designed 
to solve generalized tridiagonal eigenvalue problems but virtually, as shown 
below, it can be extended to several other rank structures. A generalization 



for tridiagonal quadratic eigenvalue problems is presented in [24| . A similar 
strategy using Newton's iteration directly applied to compute the zeros of P{z) 
is pursued in [llj. 

Modified methods for palindromic eigenproblems which are able to preserve 
their spectral symmetries have been proposed in several papers. The construc- 
tion of T-palindromic linearizations of palindromic eigenproblems is the subject 
of [Toj and Q, whereas numerical methods based on matrix iterations have 
been devised in [^^l, [3, [l^ and [2§| for computing the eigenvalues of these 
linearizations by maintaining the palindromic structure throughout the compu- 
tation. To date, however, the authors are not aware of any specific adaptation 
of the root-finding based methods to palindromic structures. 

The contribution of this paper is to fill the gap by developing a root finder 
specifically suited for T-palindromic matrix polynomials, with particular em- 
phasis on the case of large degree. T-palindromic polynomials of large even 
degree arise as truncation of Fourier series in several applications such as spec- 
tral theory, filtering problems, optimal control and multivariate discrete time 



series prediction [27 [. 



The polynomial root-finding paradigm is a flexible, powerful and quite gen- 
eral tool for solving both structured and unstructured polynomial eigenprob- 
lems. In its basic form it proceeds in four steps: 

1. The matrix polynomial is represented in some convenient polynomial basis. 

2. The transformed polynomial is linearized. 

3. The linearization is reduced in the customary Hessenberg-triangular form. 

4. A root-finding method is applied for approximating the eigenvalues of the 
(reduced) pencil. 

This scheme has some degrees of freedom concerning the choice of the polyno- 
mial basis at step 1 and the choice of the linearization at step 2 which can be 
used to exploit both structural and root properties of the matrix polynomial. 
The complexity heavily depends on the efficiency of the polynomial zero-finding 
method applied to the determinant of the pencil. Steps 2 and 3 are optional but 
can substantially improve the numerical and computational properties of the 
method. Some caution should be used at step 1 since the change of the basis 
could modify the spectral structure of the matrix polynomial. The key idea we 
propose for the implementation of step 4 is the use of the Jacobi formula [isj . 
We emphasize that, although in this paper we focus on palindromics and on a 
version of the method that is able to extract the palindromic spectral structure, 
this strategy may be used to address the most general case of an unstructured 
matrix polynomial eigenproblem, for instance by applying it to the companion 
linearization. An analysis of the application of the method to a generic matrix 
polynomial will appear elsewhere. 
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In this paper we consider the polynomial root-finding paradigm for solving T- 
palindromic eigenproblems. In particular, we address the main theoretical and 
computational issues arising at steps 1, 2 and 4 of the previous scheme applied to 
T-palindromic matrix polynomials, and also we indicate briefly how to carry out 
the reduction at step 3. The proposed approach relies upon the representation 
and manipulation of T-palindromic matrix polynomials in a different degree- 
graded polynomial basis {(f>j{y)}, namely the Dickson basis, satisfying a three- 
term recurrence relation and defined by 4>o{y) = 2, 4'i{y) = y and y(j)j{y) = 
+ (f>j^i{y) for j = 1,2,.... For the given T-palindromic polynomial 
P{z) of degree k = 2hwe determine a novel polynomial M{y) — X]^=o ■^'^j4>j{y)-: 
Mj G c'2nx2n^ 0<j<h + l, y = z + z-\ with the property that if A and 
are two distinct (i.e. A 7^ ±1) finite semi-simple eigenvalues of P{z) with 
multiplicity £, then /i = A + A^^ is a semi-simple eigenvalue for M{y) with 
multiplicity 2£. Moreover, we find that 

g{y) = det{M{y)) = [det{z-''P{z))]' ^ p{y) ■ p{y), 

where p{y) is a polynomial of degree nh at most. 

Solving the algebraic equation p{y) = is at the core of our method for T- 
palindromic eigenproblems. Our computational experience in polynomial root- 
finding indicates that the Ehrlich-Aberth method [if, Q for the simultaneous 
approximation of polynomial zeros realizes a quite good balancing between the 
quality of convergence and the cost per iteration. The main requirements for 
the effective implementation of the Ehrlich-Aberth method are both a fast, ro- 
bust and stable procedure to evaluate the Newton correction and a 
reliable criterion to stop the iteration. Concerning the first issue it is worth not- 
ing that p(y)/p'{y) = 2g{y)/g'{y) and, therefore, the computation immediately 
reduces to evaluating the Newton correction of g{y). A suitable structured lin- 
earization L{y) of M{y) can be obtained following [2J which displays a semisep- 



arable structure. In this way, in view of the celebrated Jacobi Formula|13| 
g'{y)/g{y) = trace((L(y))^^L'(?/)), the Newton correction can be evaluated by 
performing a QR factorization of L(y), say Hjj) = Q{y) ■ R{y), at low compu- 
tational cost and fulfilling the desired requirements of robustness and stability. 
Also, since || {L{y))~^ l|2=|| {R{y))~^ II2 for y ^ spec(L(y)) we obtain at no 
additional cost a reliable stop condition based on an estimate of the backward 
error given by Higham and Higham [l^ . 

If fc, the degree of the matrix polynomial, is large with respect to n, the 
size of its matrix coefficients, our approach looks appealing since with a smart 
choice of the starting points it needs 0{n^k + n^k'^) operations, whereas the 
QZ method makes use of 0{n^k^) operations. The unpleasant factor in 
our cost estimate depends on the block structure of the linearization used in 
our current implementation and can in principle be decreased by performing 
the preliminary reduction of the linearization in Hessenberg-triangular form 
as stated at step 3 of the basic scheme. The reduction can be carried out 
by a structured method exploiting the semiseparable structure of the block 
linearization to compute a rank-structured Hessenberg-triangular linearization. 
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Incorporating the structured method in our implementation would finally lead 
to a fast method that outperforms the QZ algorithm for large degrees and is 
comparable in cost for small degrees. 

The paper is organized as follows. The theoretical properties of the consid- 
ered linearizations of T-palindromic matrix polynomials expressed in the Dick- 
son basis are investigated in Section [5] and |31 The derivation of the proposed 
eigenvalue method for T-palindromic eigenproblems is established in Section |4] 
andO The complete algorithm is described in Section [51 Numerical experiments 
are presented in Section[7]to illustrate the robustness of our implementation and 
to indicate computational issues and possible improvements of our algorithm 
compared with other existing methods. Finally, conclusion and future work are 
discussed in Section [8] 

2. Theoretical preliminaries on polynomial bases linearizations 

This preparatory section recalls some basic definitions, background facts and 
notations used throughout the paper. 

For J = 0, . . . , fc let Pj G C"^", Pk ^ 0, be constant matrices and consider 
the matrix polynomial P = -P(A) = X]j=o^i'^"'- '^^^ generalized polynomial 
eigenproblem (PEP) associated to P{X) is to find an eigenvalue Aq and a corre- 
sponding nonzero eigenvector Xq satisfying 

P{\o)xo = 0. (1) 

In this paper, we will always suppose that P{X) is regular, i.e. its determinant 
does not identically vanish. 

A linearization of P{X) is defined as a pencil L{X) = XX + Y, with X,Y G 
(^knxkn^ such that there exist unimodular polynomial matrices E{X) and F{X) 
for which 

EiX)LiX)F{X) -(''[^\' ). 

Moreover, if one defines the reversal of a matrix polynomial as rev(P) := 
^ 1 the linearization is said to be strong whenever rev(L) = XY + X 

is a linearization of rev(P). 

Following the work of Mackey, Mackey, Mehl and Mehrmann [2l|, in the 
paper [T^ Higham, Mackey, Mackey and Tisseur study the two (right and 
left) ansatz vector linearization spaces: having introduced the vector A 
(1, A, . . . , A*^~^)-^, these spaces are defined as follows: 

Ci -.^ {L = XX + Y -.Bv e ^s.t.L • (A (g) /„) = V (g) P} (2) 
£,2 ■■= {L = XX + Y -.Iw e C''s.t.{A^ ®In)-L^w'^ ® P}. (3) 

It is shown in pH that almost every pencil in these spaces is a linearization, 
while in [isj two binary operations on block matrices, called column shifted 
sum and row shifted sum, are first introduced and then used to characterize the 
above defined spaces. 
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On the other hand, in ^ Amiraslani, Corless and Lancaster consider hn- 
earizations of a matrix polynomial expressed in some polynomial bases different 
than the usual monomial one. Equation (7) in Q resembles closely the defining 
equation of £2. The authors themselves stress this analogy, that suggests an 
extension of the results of [l^ to the case of different polynomial bases. Let 
{0i}i=o,...,fc be a basis for the polynomials of degree less than or equal to k. In 0] 
degree-graded bases that satisfy a three-terms recurrence relation (for instance, 
orthogonal polynomials always do so) are considered: 

X<P,{X) = a,(t>j+i{\) + l3j<f>j{X) + jj<Pj-i{X). (4) 

The aj are obviously linked to the leading-term coefficients of the . Specifi- 
cally, calling Cj such coefficients, one has that Cj = a^Cj+i. 

We wish to consider the expansion of the polynomial P{X) in this basis: 

fe 

P(A) = ^^,0,(A). (5) 

We introduce the vector 

*:= (0o(A),0i(A),...,0fe_i(A)f . 

By generalizing the linearizations studied in for each choice of $ two new 
ansatz vector linearization spaces can be defined: 

Ci := {L = XX + Y -.Bv e C^s.t.L ■ (* ® J„) = Ck-iv ® P}; (6) 
£2 := = AX + r : 3u; e C\s.t.(*^ ® I^) ■ L ^ Ck^iw'^ ® P}- (7) 

It is worth noticing that it is not strictly necessary for the new basis to be 
degree-graded, nor it is to satisfy a three-term recurrence relation. In fact, it is 
sufficient that {<j)i}i=o,...,k~i are linearly independent and have degree less than 
or equal to fc — 1, so that there exists an invertible basis change matrix B such 
that * = BK. The basis is degree-graded if and only if B is lower triangular. 
In the light of the above definitions it is immediately seen that the main 



results of 2l|, 15[ remain valid in the case of a more general polynomial basis. 
In particular the following result holds. 

Proposition 1. Let L G £1 (£2)- The following properties are equivalent: 
• L is a linearization of P 
m L is a strong linearization of P 



L is regular 



Proof. It is a corollary of Theorem 4.3 of [2l|. In fact, any L G £1 (resp., £2) 
can be written as i = Ck-iL ■ {B~^ (g) J„) (resp., L — Ck-i{B~'^ (g) /„)L) for 
some ij e £1 (resp., £2)- Therefore, L has each of the three properties above if 
and only if L has the corresponding property. □ 
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This proposition guarantees that almost every (more precisely, all but a 
closed nowhere dense set of measure zero) pencil in L\ (£2) is a strong lin- 
earization for P. For a proof, see Theorem 4.7 of [2l|. The eigenvectors of L 
are related to those of P. More precisely, (A, $ (g) a;) is an eigenpair for L if and 
only if (A, x) is an eigenpair for P . Moreover, if L is a linearization then every 
eigenvector of L is of the form ^ ® x for some eigenvector a; of P. A similar 
recovery property holds for the left ansatz vector linearizations. These proper- 
ties can be simply proved as in Theorems 3.8 and 3.14 of [2H, that demonstrate 
them for the special case $ = A. 

For the numerical treatment of palindromic generalized eigenproblems a cru- 
cial role is played by the so-called Dickson basis Q {</>i}i>o defined by 

>o(y) = 2 

<'/'i(y) = 2/ (8) 
yi > 1, v(i>3{v) = 0j+i(2/) + 0j-i(y)- 

If we consider the mapping y : = A -|- A~^ (which we will refer to as the Dickson 
transformation or the Dickson change of variable) then A^ -I- A"-' = 0j(j/) for 
j = 0,1,.... For A — e*", we obtain that (/)j(y) = 2cos(ja). From [2] by 
choosing as the ansatz vector we find a suitable strong linearization of P{X) 
represented as in ([5|): 

\ 

-/„ -7„ 

Ak-2 - Ak Ak-i J 

(9) 

In the next section we study the spectral modifications induced by the Dick- 
son change of variable that provide the basic link between palindromic matrix 
polynomials and matrix polynomials expressed in the Dickson basis. 



A + 







\ Ao 



-27„ 




3. Preservation of Jordan structure in the Dickson transformation 

Let us recall that if Aq is an eigenvalue of P(A) then the set {a;j}, = 0, . . . , £ 
is a Jordan chain of length £ -I- 1 if ajp 7^ and the following relations hold [12| : 

7 J^j — U, //t — U, . . . , t, 

^ — ' m — « 

i=o ^ ' 

where P'^'^^(Ao) denotes the fc-th derivative of P(A) evaluated at A = Aq. The 
case TO = corresponds to the definition of an eigenvector. The notion of a 
Jordan chain can be extended to any matrix function P : C — ?> C"^" whose 
determinant vanishes at Aq, as long as P(A) is analytic in a neighborhood of Aq. 
In particular, the case of Laurent polynomials is important for our investigations. 
If the principal part of a Laurent polynomial P(A) is a polynomial of degree k 
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in 1/A, then P(A) = A'^L(A) is a polynomial. The following lemma relates the 
Jordan chains of the two. The proof is a straightforward application of the 
product differentiation rule. 

Lemma 1. Let L{X) be a (Laurent) polynomial and P{X) = X''L{X) for some 
natural number k. Then the set {xj} is a Jordan chain of length £+1 for P{X) 
associated to the eigenvalue Xq ^ if and only if {xj} is a Jordan chain of 
length £ + 1 for L{X) associated to the same eigenvalue. 

Roughly speaking, Lemma [T] makes us able to switch between regular and 
Laurent polynomials without worrying about changes in eigenvalues and gen- 
eralized eigenvectors. Actually, this result can be slightly generalized with the 
next lemma, which is just an adaptation of a well-known result in 12] for the 
case where the four matrix functions that we are going to consider are polyno- 
mials. In order to prove the lemma, we recall that a vector polynomial cf>{X) 
is called a root polynomial of order £ + 1 corresponding to Aq for the matrix 
polynomial P(A) if the following conditions are satisfied: 

U(Ao)^0; ^^^^ 
I Aq is a zero of order £ + 1 for P(A)0(A). 

Obviously a root polynomial of order £+1 is defined up to an additive term of 
the form (A — Ao)^^^i'(A) for any suitable vector polynomial v{X). It is possible 
to prove that cj){X) = Ej=o(A " Ao)^^j + (A - XoY+^v{X) if and only if {cj)^} 
is a Jordan chain of length ^ -I- 1 for P(A) at A = Aq. When Aq ^ 0, thanks to 
Lemma[T]it is possible to extend the concept to Laurent polynomials: if L{X) is 
a Laurent polynomial whose singular part has degree fc as a polynomial in A~^, 
then wc say that 0(A) is a root polynomial for L{X) if it is a root polynomial 
for A'^L(A). 

Lemma 2. Let Pi (X), P2W be (Laurent) polynomials and A{X), B{X) be two 
matrix functions with P2{^) = A{X)Pi{X)B{X). Suppose that an open neighbor- 
hood Q of Xq exists such that all the considered functions are analytic in 
fl, and also suppose that both A{Xo) and B{Xq) are invertible. Then Aq is an 
eigenvalue for Pi if and only if it is an eigenvalue for P2, and {y^} is a Jordan 
chain of length £ + 1 for P2 at Aq if and only if {zi} is a Jordan chain of length 

£+1 for Pi at Xq, where Zi = J2]=o ^—^Vi-j- 

Proof. If Pi{X) and -P2(A) are classical polynomials then the thesis follows 
as in the proof of Proposition 1.11 in [l2| after having represented ^(A) and 
B{X) by their Taylor series expansions. To deal with the Laurent case, let 
a and /3 be the minimal integers such that Qi{X) := A"Fi(A) and Q2{^) ■= 
X^P2{X) are classical polynomials. Just follow the previous proof for Q2{X) = 
A'5-"A(A)Qi(A)B(A) and apply Lemma [HD 

We are now in the position to prove a result for the Dickson change of variable 
y = X + l/X. The following proposition shows that the number of Jordan chains 
and their length at some eigenvalue yo (for the sake of brevity, we shall use the 
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expression Jordan structure at yo) is related to the Jordan structures at Aq and 
A-. 

Lemma 3. Let y{\) = A + A^^ and let M{y) be a polynomial in y, so that 
7V(A) :~ AI{y{\)) is a Laurent polynomial in A. Let first j/o = Ao + A(J^, 
Ao 7^ ±1, be a finite eigenvalue of M{y). Then the Jordan structure of M{y) 
at yo is equal to the Jordan structure of N{X) at either Ao or I/Aq. // on the 
contrary Aq = ±1, then there is a Jordan chain of length i at M(±2) if and 
only if there is a Jordan chain of length 2£ at A'^(±l). 

Proof. It is obvious that yo £ C is an eigenvalue for M{y) if and only if 
both Ao and A,^^ are eigenvalues of N{X). Let M(y) = E( y)D(y)F(y), where 
D(y) — diag{di (y ),... ,dn{y)) is the Smith form ([l3|,[2a|) of M{y). Define 
^(A) := E{ylx)),DiX) D{y{X)),F{X) := F(y(A)). If a,/3,7 are such that 
N(X) = A"+'5+TiV(A), E{X) = A"£;(A), D{X) = A''l)(A) and ^(A) = A''F(A) 
are polynomials in A, then we have the relation N{X) ~ E{X)D(X)F{X); how- 
ever, in general D{X) is not the Smith form of N{X). Nevertheless, it has the 
form diag(A'^i(ii(A), . . . , A*''"d„(A)) where fci > fc2 > • • • > A:„ and di{X) = 
yicg{di)^_^^y^-^jy jjj Other words, the c?i(A)-s are palindromic polynomials with 
no roots at and such that di{X) divides di+i (A) for i — 1, . . . ,n~ 1. Moreover, 
yo is a zero of multiplicity n for di{y) if and only if both Ao and Aq ^ are zeros 
of multiplicity n for di{X). To reduce D{X) into a Smith form, we proceed by 
steps working on 2 x 2 principal submatrices. 

In each step, we consider the submatrix '^'^^^ ^b!; /^^. h with i < j. If 

y A dj (A) J 

a < /3, then do nothing; if a > /3, premultiply the submatrix by ( _f,^(A) i-fc(A) ) 
and postmultiply it by (^^^^^^ a"'^^ ) ' ^^^'''^ lW = '^jW/'^iW while a(A) and 
6(A) are such that a(A)A"c?i(A) + b{X)X'^ dj{X) = X^di(X); the existence of two 
such polynomials is guaranteed by Bezout's lemma, since X^di{X) is the great- 
est common divisor of X°'di{X) and X^dj{X). It is easy to check that both 
matrices are unimodular, and that the result of the matrix multiplications is 

d,(A) ^^9^^^^. By subsequent applications of this algorithm we thus con- 
clude that the Smith form of D{X) is D{X) = diag(A'="(ii(A), . . . , A'=iJ„(A)). 

It follows that the ith invariant polynomial of M{y) has a root of multiplic- 
ity n.i at yo if and only if the ith invariant polynomial of D{X) has a root of 
multiplicity at Ao ±1 and a root of multiplicity Ui at I/Xq. From Lemma 
[21 the Jordan structures of iV(A) are equal to those of D{X). The thesis follows 
from the properties of the Smith form and from Lemma [TJ 

Mutatis mutandis, a similar argument can be used to analyze the case of 
A — ±1: notice in fact that {y ± 2)'^ is a factor of the ith invariant polynomial 
of M{y) if and only if (A ± 1)^'^ is a factor of the ith invariant polynomial of 
D{X). □ 
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4. Application to palindromic polynomials 

We will now specialize our analysis to the case of a matrix polynomial with 
palindromic structure. 

Remark 1. In this section, we will only treat the case of even degree palin- 
dromic matrix polynomials. Notice in fact that an odd degree palindromic may 
always be transformed to an even degree palindromic, either by squaring the 
variable (A = /x^) or by multiplication by (A + l)/„. Potentially, both actions 
may introduce problems: squaring the variable adds an additional symmetry 
{/i, — /i} to the spectrum while multiplying by A + 1 increases by n the multi- 
plicity of —1 as an eigenvalue. 

However, the first issue may be solved, after passing to Laurent form, by the 
use of the change of variable z = (/i + /i^^)^. See also Remark [3l 

Regarding the latter issue, since one knows that he is adding n times — 1 
there is no need to compute it: n of the (n + l)k starting points of the Ehrlich- 
Aberth iteration shall be set equal to —2, and there they remain with no further 
corrections. The shortcoming is that the Jordan structure at A = — 1 changes. 

Let P(A) = Y^'^j'Lo be a polynomial of even degree. By Lemma [1] 

switching to the Laurent form is not harmful for finite nonzero eigenvalues 
and the corresponding (generalized) eigenvectors; we can therefore consider its 
Laurent counterpart 

fc 

P(A) := ^ A,X^. (11) 

Three different kinds of palindromic structure can be defined. We say that 
the Laurent polynomial is purely palindromic (resp., ★-palindromic, ★ G {T, H}) 
if the following relations hold between its matrix coefficients: 

Purely palindromic: Aj ~ A^j] 
★-palindromic: Aj = A*_j. 

It is well-known that the palindromic structure induces certain symmetries of 
eigenvalues and eigenvectors: in particular if Aq is an eigenvalue, a; is a right 
eigenvector and is a left eigenvector, then, denoting complex conjugation 
with the operator (•)* 

if P is purely palindromic, P[j^)x = 0, P[j^) = 0; 
< if P is T-palindromic, P{^)z = 0, x^P{^) = 0; 
if P is i/-palindromic, P{^)z* = 0, x" P{^) = 0. 

In this paper we are primarily interested in the design of an efficient solver 
for T— palindromic eigenproblems. A numerical method will be presented in 
Subsection 14.21 The proposed approach can however be described very easily 
with purely palindromic polynomials. Thus we first consider this case for the 
sake of clarity. 
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Purely palindromic polynomials 

The most obvious way to deal with this kind of pahndromicity is via intro- 
duction of the change of variable y = A + A~ ^ , in order to halve the degree of 
the polynomial. More explicitly, one can define Q{y) :~ P{X{y)); clearly, the 
purely palindromic structure of P{X) guarantees that Q{y) is itself a polynomial 
in the new variable y. The next proposition is a simple application of Lemmas 
[T] and [31 and it relates eigenvectors and Jordan chains of the two polynomials: 

Proposition 2. When Aq ± 1, the Jordan structure of Q{y) at the eigenvalue 
yo = -^0 + Aq"^ is equal to the Jordan structure of P{\) at either Aq or Ag ^- If 
Ao = ±1, Q{y) has a Jordan chain of length i at yo — ±2 if and only if P(X) 
has a Jordan chain of length 2£ at Ao = ±1. 

In particular, the eigenvectors of Q{y) at yo ai'e exactly the same of the 
eigenvectors of P(A) at Aq (or equivalently at A,^^, since they are the same). 

Albeit very attractive, from a numerical point of view this trick is not very 
suitable as soon as one considers a high degree polynomial. In fact, the matrix 
coefficients of Q(y) need to be computed as linear combinations of the ones 
of f (A). Since the powers of a binomial are involved, the coefficients of these 
linear combinations would exponentially grow with the polynomial degree. To 
circumvent this difficulty, we shall make use of the Dickson polynomials ^ . The 
polynomial Q{y) is readily expressed in terms of the 4>j{y)s since in the Dickson 
basis the coefficients are just the old ones and therefore no computation at all 
is needed, namely, 

4 * 

Q(y)-^'/'o + E^^'^j(y)- (12) 

The associated linearization ([9|) has several computational advantages with re- 
spect to other customary linearizations of P{X). Its size is nk versus 2nfc, the 
spectral symmetries are preserved and, moreover, the linearization displays a 
semiseparable structure. More precisely, it is of the form Dq + Diy where Di 
is identity plus low rank while Dq is Hermitian plus low rank. This kind of 
structure is preserved under the QZ algorithm and it may be exploited for the 
design of an efficient and numerically robust root-finder applied to the algebraic 
equation det Q{y) = 0. 

4-2. T-palindromic polynomials 

Consider now a T-palindromic polynomial of even degree 2k. We will sup- 
pose once more that neither nor oo are eigenvalues, so that we can divide 
by A'^ and consider the Laurent form P{X), which is a T-palindromic Laurent 
polynomial of degree k both in A and in A~^. Since the symmetry X -(r^ X~^ is 
still present in the spectrum, we expect that the Dickson basis may still play a 
role. However, unlike the purely palindromic case, it is not possible to directly 
express a T-palindromic polynomial as a polynomial in the variable y. In fact. 
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splitting P{X) as the sum of its symmetric part and its skew-symmetric part we 
obtain that 



P{X)=Ao + J2 



A, + Ai 



A. 



(13) 



If we introduce the new variables y := A + A~^ and w :— X — X~^ , then P(A) 
can be expressed as a bivariate polynomial in w and y which is always linear in 
ui, that is, 

Q{y, w) = P{X{y, w)) B{y) + wC{y). 



The property follows from (1131) by substituting 

. 1 + (-iy+^ 
A^ +A-J =0j(2;), A^ -A-^ =- ' ^ ' 



j>i- 



Notice moreover that B{y) is a symmetric polynomial (that is to say, every 
matrix coefficient is symmetric), C{y) is skew-symmetric, and the operation of 
transposition corresponds to changing the sign of w, that is, 

g^(2/, w) = P^{X(y, w)) - B{y) - wC{y). 

In principle one may think of treating Qiv, w) with available techniques for the 
bivariate eigenvalue problem (see e.g. [l6| and references therein), but actually 
y and w are not independent. They are related by the trigonometric dispersion 
relation w"^ ~ y"^ — 4. This suggests that it is possible to obtain a univariate 
polynomial by doubling the dimensions of the matrix coefficients. Let us define 



Miy) 



B{y) w^C{y) 
C{y) B{y) 



Then M{y) is a polynomial in y of degree fc -I- 1 at most. Moreover, it has 
the following property: if Aq and A^"'^ are two distinct (i.e. Aq 7^ ±1) finite 
semisimple eigenvalues of P(A) with multiplicity m, then yo = ^^o + Aq"^ is a 
semisimple eigenvalue for M{y) with multiplicity 2to. To see this, notice first 
that 



Miy) = diag(y^/„, i/v^/„) ( ^^y^^ ""^l^y^ 



diag(l/Vu'-^n, Vwin) 



and 

/ B{y) wC{y) 
\ wC{y) B{y) 

Hence, we find that 
M{y) = E{w) 



Q{y,w) 




-^71 

^71 I71 







Q{y,w) 






Q^iy,w) 
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Since, as long as E{w) is defined (that is to say w ^ 0, oo or A 7^ 0, ±1, cjo), 
det(-B(w)) = 1 then 

det(Af(2/)) = [dct(Q(y,w))]2, V (y, u;) € C x C. (14) 

Therefore, Aq has algebraic multiplicity m for P(A) if and only if yo has algebraic 
multiplicity 2m for M{y). This gives the factorization 

det(M(y))=p(y).p(2/), (15) 

for a suitable polynomial p{y) having the zero yo of multiplicity m. Concern- 
ing eigenvectors, if Aq is semisimple, then let Xj (resp. Zj), j — 1, . . . ,m be 
the eigenvectors for P(A) (resp. P'^(A)) corresponding to Aq: it can be easily 
checked that {{wqxJ , xj)'^ , {—wozj, zj)'^}, where wq = \o + Xq^ , are two lin- 
early independent eigenvectors for M{y) corresponding to yQ. Thus, geometric 
multiplicity is also 2m. Indeed, something more can be said in the more general 
case of Jordan chains. 

Proposition 3. Let yo = Aq + X^^ be an eigenvalue of M{y) so that Ao and 
Aq^ are eigenvalues for P{X). If Xq ^ 0, ±l,oo then the Jordan structure of 
M(y) at yo is equal to the union of the Jordan structures of P{X) at Aq and at 

Proof. Since P(A) is T-palindromic, it is clear that the Jordan structure of 
R(X) - ( 



Aq ^ . Define 



at either Ao or Aq ^ is the union of the Jordan structures of P{X) at Ao and at 

iV(A) M(y(A)) E{w{X))R{X)E'^ {w{X)). 

The matrix function E{w), defined in the previous page, is analytic everywhere 
in the w complex plane but on a branch semifine passing through the origin. 
Since by hypothesis 7^ 0, the branch cut can be always chosen in such a way 
that E{w) is analytic in a neighborhood of wo = Xq — Xq^ , and thus E{w{X)) is 
analytic in a neighborhood of Ao- Then we can apply Lemma [5] to conclude that 
the Jordan structures of Af(A) and i?(A) are the same. Application of Lemma 
[3] completes the proof. □ 

Remark 2. Another remarkable property of M{y) is that its coefficients are all 
skew-Hamiltonian, that is to say they can be written as JK where J = C^j J) 
and K is some skew-symmetric matrix. This link between T-palindromic and 
skew-Hamiltonian polynomials is interesting because it may shed more light on 
the relation between several polynomial structures. It is known that one can 
easily transform a palindromic polynomial to an even polynomial by a Cayley 
transformation, and then to a Hermitian polynomial via a multiplication by i (if 
one started from a real polynomial) or to a symmetric polynomial by squaring 
the matrix coefficients. On the other hand, Hamiltonian polynomials can lead to 
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skew-Hamiltonian polynomials by squaring each coefficients, and multiplication 
by J sends a skew-Hamiltonian polynomial to a skew-symmetric polynomial. 
The Dickson change of variable, followed by doubling the dimension, is able to 
map T-palindromic polynomials of even degree to a subset of skew-Hamiltonian 
polynomials. Unlike some of the other mentioned maps, this is not a bijection 
between two classes of structured polynomials, because what is obtained is ac- 
tually a subset of skew-Hamiltonian polynomials. In fact, since the north-west 
and south-east coefhcients of M{y) are the coefficients of B{y) they must be 
symmetric and there is a relation between the north-east and south-west coef- 
ficients of M{y). However, a deeper investigation on this subject is needed in 
the future. 

Remark 3. Notice that a similar technique can be applied to even/odd matrix 
polynomials, that is polynomials whose coefficients alternate between symmet- 
ric and skew-symmetric matrices. In this case, on can apply the transformation 
z = and use algebraic manipulations, akin to the ones described for the 
T-palindromic case, in order to build a new polynomial in z with double dimen- 
sions. 

In the case of an odd-degree T-palindromic polynomial, the substitutions 
X — fj,^ and y — fi + fi~^ lead to an M{y) such that ( j q) ' ^iv) '^^ odd. 
Therefore, one may apply z = y'^ and build a third polynomial in order to 
extract the additional structure {/i, —fi}. 

Equation and (fT5|) enable the computation of the eigenvalues of P(A) 
to be reduced to solving algebraic equations. From Proposition |3] it follows that 
possible discrepancies in the Jordan structures can be expected for j/g — i2 and 
yo = oo corresponding to Aq = ±1 and Aq = 0, oo, respectively. 

When Ao — ±1 not only the proof we gave is not valid (because, since wq = 
is a branch point, there is no neighborhood of analyticity of the matrix function 
E), but in fact the proposition itself does not hold. As a counterexample, let 
a 7^ ±^ and consider the polynomial 



\-2 + X~^ aX-aX-^ 
~aX + aX^^ A-l-A"^ 

We have that {(1,0) ,(0,a) } is a Jordan chain for P{X) at A = 1. The 



P(A) 

corresponding M{y) is 
M{y) = 



/ y - 2 ay^ _ 4a \ 

y Aa- ay'^ 

a y~2 

V -a y 1 



which has a semisimple eigenvalue at y = 2 with the corresponding eigenvectors 
(0,0, 1,0)"^ and (2,0,0,af . 

If the leading coefficient of P{X) is not symmetric, then M (y) has 2ti extra 
infinite eigenvalues, where n is the dimension of the matrix coefficients of P{X). 
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These eigenvalues are defective since their geometric multiplicity is only n + 
dimkerCfc_i, where Ck-i is the leading coefficient of C{y). 

For the numerical approximation of the roots of p{y) we can exploit again 
the properties of the Dickson basis to compute the matrix coefficients of M{y) = 
Sj=o ^^j4>j{v)- The code below computes the matrices Mj e C^"^^", < j < 
fc — 1, given in input the coefficients Aj of P(A), < j < k, defined as in pT|). 

function Dickson_transform 

Input: Aa,...,Ak G C"^" 
Output: Mo, Mk+i e C2"^2« 

Bf) = Ao/2; Co — 0„; 

for j = 1, . . . , k 

B, = {A, + Aj)/2; C, = {A, - Aj)/2; 

end 

'S'o = 0„, 5*1 — 0„; 
for j = fc, . . . , 1 

'^modo,2) ~ '^mod(i,2) + 
^J-i ^ '^mod(j,2)' 

end 

Co = Co/2] Ck = Cfe+i = 0„; Co = C2 
Ci = Ci + C3; C2 — 2Ci + C4; 
for j — 4, . . . , k 

Cj-i = Cj-3 + Cj+i\ 

end 

Cfc = Ck-2\ Ck+1 = Cfe_i; 
for j = 1: fc + 2 

Cj-i = Cj-i - 2Cj_i; 

^j-i = [-Sj-i,Cj_i;Cj_i,Bj_i]; 

end 

Remark 4. The coefficients of C{y) are linear combinations of Aj —Aj. As can 
be seen by the above algorithm, the coefficents of such combinations expressed in 
the Dickson basis remain bounded, the upper bound being 1/2. An analogous 
result, with upper bound 1, holds for w'^C{y). This is in contrast with the 
exponential growth that would have been seen in the purely palindromic case 
if one had directly applied the Dickson transformation without the use of the 
Dickson basis. 

The arithmetic cost is 0{n?k) operations. Once the coefficients Mj are 
determined, a linearization of M{y) of the form © can be constructed. The 
properties of this linearization are investigated in the next section in order to 
devise a fast and numerically robust method to evaluate the Newton correction 
of p{y) defined by . 
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5. Computing the Newton correction 



Our aim in this section is to derive a fast, robust and stable method for com- 
puting the Newton correction p{y)/p'{y) = 2det(Af(y))/(det(Af(i/)))', where 
p{y) and M(y) are related by (jlSp . given a structured linearization L{y) = 
yE + F, with E,F e C'^^ik+i)x2n{k+i) ^ ^^ic form ©, namely, 



E : 



i2n 



Mfe+1 / 



and 



/ 



F : 





-l2n 



-l2n 




Mk-2 





A4_i - 



Our approach relies upon the celebrated Formula of Jaco6z[13j 

(det(L(2;)))' = det(L(y))trace(i-i(2/)i'(2;)) - det(i(y))trace(L-i(y)£;) 

which reduces the evaluation of det (M (y ) ) / (det (M(i/) ) )' = det (L (y) ) / (det (L (y) ) )' 
to computing the trace of L~^{y) ■ E. In the sequel we describe a method for 
finding the block entries and, a fortiori, the trace of the inverse of L{y) from 
the LQ factorization of the matrix. Then we slightly modify the computation 
to take into account the contribution due to the matrix E. It will be clear 
from what follows that this method is general and can be applied, with only 
trivial modifications, to any kind of unstructured matrix polynomial, simply by 
considering for instance the standard companion linearization instead of 
We denote as 0{d,ij}) the 2x2 unitary Givens rotation given by 



\9\' + \^j\' = l. 



Let L{y) ^ L ■ Q he the (block) LQ factorization of L{y) obtained by means of 
Givens rotations so that 



(16) 



It can be easily checked that the lower triangular factor L has the following 
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structure 



L 



( aihn 

Plhn 
llhn 



a2hn 

Ml 



Mk-2 Mk-1 Mk J 



where aj 7^ 0, 1 < < fc. If Mk and, therefore, L{y) is invertible then the 
LQ factorization can be used to find a condensed representation of the inverse 
of L{y). Observe that L^^{y) = ■ L~^. In order to take into account the 
occurrence of the matrix E in the Jacobi formula let us introduce the matrix 
Mk+i = Mj^^ ■ Mk+i- Then we have the following 



Proposition 4. There exist matrices AIi 



L-\v)E. 



/ Ml ^_iM2 
01 M2 



\ 



. Vi ■ ■ • tpkMk+i \ 



OkMk^ 



where the blank entries are not specified. 

Proof. The proof basically follows by applying the (block) Schur decompo- 
sition ([T6l) of to the block lower triangular factor L~^E. To show it more 
formally we can proceed by induction. Let us assume that the the j— th block 
row of Gj ■ ■ ■ GkL^^E can be represented as 



[ 



ipj ■ ■ ■ il^kM) 



k+l 



where Mj is the diagonal entry and the value of the entries in the strictly lower 
triangular part - denoted by is not essential. Then, by applying Gj-i on the 
left of the matrix we find that the (j — 1) — th block row looks like 



[ 



M. 



ipj-iMj 



^jjj-i ■ ■ ■ -ipkMk+i 



whereas the diagonal entry in position j becomes 9j-iMj. □ 

This result says that the block diagonal entries of L~^{y) can be determined 
from the entries in its first (block) row. The computation of this row is equivalent 
to the solution of the linear system 



{hn, 02„, . . . , 02„) = (Xi, . . . , Xk+i) ■ L{y) 



or, equivalently, 



{l2m 02n, . . 



,02 



{Xi, . . . , Xfc+i) • L. 
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In the view of the structure of this reduces to 

fe- 1 k \ 

i=i i=i / 

Let D e C'^n(k+i)x2n(k+i) ^^Jq^Jj diagonal matrix defined by 

fc-i fe 

r> = diag(l, V"!, • • • , ]J i^j, Yi ^j) ® -^2n- 

Using the matrix D to balance the coefficient matrix yields 

(^1, ^2, . . . , ^fe, 1) /2n = (^1, . • • , ^fe+i) D-^-D-L-D-\ 

Observe that 

(Xi, . . . ,lfe+i) = (Xi, . . . , Xfc+i) = (Mi, . . . , Mk-uMk^) , 
and, therefore, the solution of 

(0i,02,...,^fe,l)®/2„= (Xi,...,Xfe+i)L, DLD-^=L, 

gives the desired unknown matrices Mi, . . . , M^. To achieve some computational 
savings we rewrite the system as 

{91,92, ...,9k,l)®Mk = (^1, . . . ,^fe+i) L 
and thus we arrive at the following relation 

det(M(2/))7det(M(|/)) = trace(M-i(Xi + ^1X2 + . . . + ^fc-iXfc + hMu+i)), 

which is used to compute the reciprocal of the Newton correction. The function 
trace below implements our resulting algorithm at the cost of 0{nPk + n^) 
operations. 

function trace 

Input: Mo, . . .,Mk+i G C2"x2", A e C, (det(M(A)) ^ 0) 
Output: the value of 77 = p'{X}/p{X} 

Mk-i - Mk-i - Mk+i: Mk = Mfe + AMfc+i; 

a = A ones(l, fc + 1); /3 = —ones{l, k + 1); 

7 = zeros{l, k); x = —ones{l, k + 1); xi = ~2; 

for j = l,...,fc 

^ = K; Xj]; = planerot{v); q{j, : ) = ^(1, : ); cj = q{j, 1); 

ctj = ajGi.i + XjQ2,i; P = lijQi,i + 01J+1Q2X: 

= PjG^a + 0^+1^/2,2; = 13; 7j = /3j+iG2,i; Pj+i = Pj+iG2,2\ 
M = gi,iMj_i + g2,iMj; Mj = gi,2Mj_i + g2,2Mj; Mj_i = M; 
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end 

for j = 1, . . . , fc — 1 

/?, =/3jq(j,2); 

end 

ioi j — I, . . . ,k — 2 

7j =7j9(i,2)q(j + l,2); 

end 

ioi j ^ k, . . . ,1 

s^sq{j,2); Mj_i = sMj_i; 

end 

Xk = {ckMk - Mk-i)/ak\ Xk-i - {ck-iMk - Mfe-2 - h-M/ak-i] 
for j = fc — 2. . . . , 1 

Xj = {cjMk - Mj-i - l3jXj+i - 7jXj+2)/aj; 

end 

M = Xi- 

for J = 1 , . . . , fc — 1 

M = M + CjXj+i; 

end 

M = M + CfcMfe+i; M Mk\M; ry = trace(M); 

6. The Ehrlich-Aberth algorithm for T-palindromic eigenproblems 

A simple tool for the simultaneous approximation of all the eigenvalues of a 
polynomial is the Ehrlich-Aberth method. Bini and Fiorentino [4] showed that 
a careful implementation of the method yields an efhcient and robust polyno- 
mial root finder. The software package MPSolve documented in 0] is designed 
to successfully compute approximations of polynomial zeros at any specified 
accuracy using a multi-precision arithmetic environment. 

A root finder for T-palindromic eigenproblems can be based on the Ehrlich- 
Aberth method applied for the solution of the algebraic equation p{y) ~ 0, 
where p{y) is related with M{y) by ([T5)) and M{y) is generated by the function 
Dickson_transform applied to the input coefficients Aj € C"^" of the T- 
palindromic matrix polynomial -P(A) of degree 2k given as in (IT5|) . The method 
simultaneously approximates all the zeros of the polynomial p{y): given a vector 
G C^, N = nk, of initial approximations to the zeros of p{y), the Ehrlich- 
Aberth iteration generates a sequence {z^'^)}, fc > 0, which locally converges to 
the A^— tuple of the roots of p{y), according to the equation 

The convergence is superlinear for simple roots and linear for multiple roots. 
In practice, the Ehrlich-Aberth method exhibits quite good global convergence 
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properties, even though no theoretical results are known in this regard. The 
main requirements for an efficient implementation of the method are: 

1. a rule for choosing the initial approximations; 

2. a fast, numerically robust and stable method to compute the Newton 
correction p{z)/p'{z); 

3. a reliable stopping criterion. 

Concerning the first issue it is commonly advocated [s!] that for scalar poly- 
nomials the convergence benefits from the choice of equally spaced points lying 
on some circles around the origin in the complex plane. In the case of matrix 
polynomials where the eigenvalues are often widely varying in magnitude this 
choice can not be optimal. A better strategy using the initial guesses lying 
on certain ellipses around the origin in the complex plane is employed in our 
method. The second task can be accomplished by means of the function trace 
in the previous section. With respect to the third issue, it is worth observing 
that the QL-based method pursued for the trace computation also provides an 
estimate on the backward error for the generalized eigenvalue problem. From a 



result in it follows that if y is not an eigenvalue of L{y) then 
rKy) = 1/(11 {yE + F)-^ Il2(l + |y|)) 

gives an appropriate measure of the backward error for the approximate eigen- 
value y. Since for yE + F = L ■ Q we have that 

II (yE + F)-' h=\\ |h>|| M^' \\2> iV2^)-' II M^:' lico . 

In our implementation we consider the quantity 

fj{y)^V27i/i\\M-' l|oo(l + |y|)) 

as an error measure. If ri{y) is smaller than a fixed tolerance then y is taken 
as an approximate eigenvalue and the corresponding iteration is stopped. The 
resulting Ehrlich-Aberth algorithm for approximating finite eigenvalues of M{y) 
and hence obtaining the corresponding eigenvalues of P{X) is described below. 
In the next section we present results of numerical experiment assessing the 
robustness of the proposed approach. 

function palindromic_aberth_zeros 

Input: Ao, . . . , Afe G C"^", tol E R, maxit £ N, initial guesses zi, . . . ,zn 
Output: approximations Ci, . . . , (2N, N — nk, of the zeros of P{X) = X]i=-fc ^j'^' 
[Mo, . . . , Mk+i] =Dickson_transform(Ao, . . . , A^) 

TV = nfc; c = ones(N, 1); 
nn = 0; 

for i = 1, . . . , maxit 

for j = 1,...,N 

if m) 

z =trace(Mo, . . . , Mk+i, zj); z = 2/z; 
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h = sum{l./{z{l:3-l)-z{j)))- 

h = h + sum{l./{z{j + 1 : N)- z{j))); 

h — z/{l + h z); Zj = Zj — h] 

if {'fi{zj) < tol or \h\ < tol\zj\) 

c{j) = 0; nn = nn + 1; 

end; 

end 

end 

if {nn ^ N) 
break 

end 

end 

forj^l,...,N 

r — roots{[l, ~Zj, 1]); 

end 

The total cost of the algorithm is therefore 0{t{n'^k + n^)) operations, where 
t is the total number of times that the function trace is called. Numerical exper- 
iments presented in the next section show that t heavily depends on the choice 
of the starting points. With a smart choice, t is of order 0{nk), which gives a 
total computational cost of 0{n'^k + n^k'^). Since the cost of our method grows 
as n'^ but is only quadratic in fc, where customary QZ-like methods use 0{n^k^) 
operations, an Ehrlich-Aberth approach looks particularly suitable when the 
matrix polynomial has a high degree and small coefficients so that k'^ /n is large. 

It is worth noticing that the case of large n can still be treated by means of 
an Ehrlich-Abcrth method in 0{n^k^) operations. The basic observation is that 
the factor n"* comes from the block structure of the linearization involved in the 
computation of the trace. A reduction of the cost can therefore be achieved by 
a different strategy where the linearization is initially converted into (scalar) 
triangular-Hessenberg form: say, N{y) = Ry + H where R is (scalar) triangular 
and H is (scalar) Hessenberg. The task can virtually be performed by any 
extension of the fast structured methods for the Hessenberg reduction proposed 
in . These methods preserve the rank structure which can therefore be 

exploited also in the triangular-Hessenberg linearization. Once the matrices R 
and H have been determined then the computation of tv{N{y)) can be performed 
by the following algorithm which has a cost of 0{n^k^) = 0{N'^) operations: 

• Perform a RQ decomposition of the Hessenberg matrix N{y), obtaining a 
unitary matrix Q represented as product oiO{N) Givens transformations 
(Schur decomposition) and a triangular matrix U. 

• Compute the last row of N{y)^^R by solving w'^N{y) — ej^ and then 
computing w'^ : — R. 

• Recover the diagonal entries of N{y)^^R from the entries of w and the 
elements of the Schur decomposition of Q. 
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This alternative road leads to an algorithm of total cost 0{n^k^) operations. 
An efficient implementation exploiting the rank structures of the matrices in- 
volved will be presented elsewhere. 



7. Numerical Experiments 

The function palindromic_aberth_zeros for computing the roots of a T- 
palindromic matrix polynomial P( A) = X] j=-fe i given its coefficients A^j = 
Aj e C"''", < j < fc, has been implemented in Matlalfl and then used for 
the computation of the zeros of polynomials of both small and high degree. The 
tolerance is fixed at tol = l.e — 13 and for the maximum number of iterations 
we set maxit = 2nk. 

Extensive numerical experiments have been performed to illustrate some 
basic issues concerned with the efficiency and the accuracy of a practical imple- 
mentation of our method. 



7.1. Efficiency of root-finding 

An accurate and efficient root-finder is essential to the success of our al- 
gorithm. In practice, the cost of each iteration is strongly dependent on the 
amount of early convergence (for the sake of brevity, in the following we will 
refer to this phenomenon using the word deflation) occurring for a given prob- 
lem. In other words, a critical point to assess the efficiency of the novel method 
is the evaluation of the total number t of calls of the function trace, and of its 
dependence on the total number N := nk of the eigenvalues. When the Ehrlich- 
Aberth method is used to approximate scalar polynomials roots, experiments 
show that t depends on the choice of the starting points. If there is not any a 
priori knowledge about the location of the roots, empirical evidence Q shows 
that choosing starting points distributed on some circles around the origin leads 
to acceptable performances and/or quite regular convergence patterns. 

The class H„_fe of T-palindromic polynomials have been used to verify if 
these properties still hold in the matrix case. The polynomials are constructed 
according to the following rules: 

Ao = 0„; Aj = /„ + e„e^, A^j ^ Aj , l<j< k. 



From 



3 



Xk _ 1 A'^+i -t- 1 
A^ ' 



we find that most of the eigenvalues lie on the unit circle and for k even A = — 1 
is a double root of h{X). 



^Matlab is a registered trademark of The Mathworks, Inc.. 
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Figure 1 describes the convergence history for our root finder applied to 
-^5,20 with starting values equally spaced on the circle centered in the origin 
with radius 4. The curves represented are generated by plotting the sequences 
{^j'^^l, 1 < fc < maxit, for j = 1, . . . , A^. The convergence is quite regular and 
very similar to that exhibited in the scalar polynomial case and theoretically 
predicted for simultaneous iterations based on Newton- like methods fv^. 




Fig. 1. History of the convergence for the H problem with n — 5 and k — 20 



With this choice of starting points we have observed that the number of 
global iterations is typically of order of N but there are not enough early defla- 
tions, that is, iterations that are prematurely stopped due to early convergence. 
In order to increase the cost savings due to premature deflation in our program 
we have employed a slightly refined strategy. Since the method does not ap- 
proximate directly the eigenvalues but their Dickson transform yi = Xi + X^^, 
we have chosen starting points on the Dickson transform of the circles \z\ = p, 

that is points lying on ellipses ^^j^^^^ -I- (p™|/p)2 = 1- More precisely, this is 
the algorithm we used to pick the starting points: 



Input: Number N of eigenvalues to approximate and parameters a G N and b E N 
Output: Starting points Zk, k = 1, . . . , N 

= 2tt/N- 

(/)=randn; 

for j = 1,...,7V 

jj=mod(j,a); 
P = 1 - ji/b; 
a = p+l/p\ 
P=l/p-p- 

Zj — a cos{j * 6 + 4>) + (3 sin(j * 9 + (f)) 
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end 



The integer a determines the number of ellipses whereas b is used to tune 
the lenghts a and /3, defined as above, of their semiaxes. We expect that a good 
choice for the parameters a and b depends on the ratio k/n: when k ^ n we 
expect many eigenvalues to lie on or near to the unit circle, while when n ^ k 
we expect a situation more similar to the eigenvalues of a random matrix, with 
no particular orientation towards unimodularity. We therefore expect that a 
small ratio a/b works well in the former case while on the contrary in the latter 
case a ~ 6 should be a better choice. Moreover, we expect that as nk grows it 
is helpful to increase the total number a of ellipses as well. 

We show here some of the results on random T-palindromic polynomials. 
Figure 2 refers to an experiment on small-dimensional, high-degree polynomials: 
the value of n has been set to 5 while k was variable. The average number of 
t over a set of 1000 random polynomials for each value oi N ~ nk is shown 
on the graph. The parameters satisfy a £ {2,3} and b G {8,64} and they 
are determined by a = 1 + 2"^ and b = 8'^^^, where the integer c is defined as 
c = lo 

§320 ^- The graph shows a linear growth of t with respect to iV — nk. 
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Fig. 2. 

Figure 3 refers to an experiment where on the contrary the case of small k is 
explored. We have considered here k = 2 and let n vary and we show the results 
for t plotted against nk for several choices of a and b. The choice labelled as 
'step function' is for a — {6, 11} and b — {6, 12} generated by a = 1 -I- 5 2"^ and 
6 = 62"^. Once again the experiments suggest that when the starting points 
are conveniently chosen t < aN for some constant a and any N in the specified 
range, and, moreover, the bound still holds for different reasonable choices of 
the parameters a and b. The experimentation with random polynomials gives 
a ~ 8 as an estimate for the constant. 



n=5 
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Fig. 3. 

In conclusion, the algorithm can greatly benefit from a smart strategy for 
the selection of the starting points by increasing the number of early deflations. 
The experiments show that as long as the starting points are suitably chosen 
the value of t is proportional to N — nk. 

7.2. Accuracy of root-finding 

The other important aspect of our solver based on polynomial root-finding 
concerns the accuracy of computed approximations. In our experience the 
method competes very well in accuracy with the customary QZ-algorithm. The 
accuracy of the computed non-exceptional roots for the random polynomials 
was always comparable with the accuracy of the approximation obtained with 
the QZ method. The results of other numerical experiments confirm the ro- 
bustness of the novel method. Figure 4 illustrates the computed eigenvalues for 
the problem Figure 5 also reports the plot of the absolute error vector 

cihs{XEA~ ^) and abs{\Qz — X), where A is the vector formed by the eigenvalues 
computed in high precision arithmetic by MathematiccH while A^^i and Xqz 
are, respectively, the vectors formed from the eigenvalues returned by our rou- 
tine paIindromic_aberth_zeros and suitably sorted by the internal function 
polyeig. 



^Mathematica is a registered trademark of Wolfram Research, Inc. 



24 



1.5 




-1.5 -1 -0.5 0.5 1 

Eigenvalues for the H problem with ri-5 and k=40 



Fig. 4. 




Fig. 5. 

The numerical results put in evidence the following important aspects: 

1. Poor approximations for the exact eigenvalue A = — 1 are in accordance 
with the theoretical predictions: in fact the reverse transformation from 
y = X + X^^ toA = ^(y± — 4) is known to be ill-conditioned near 
y = ±2 (or A = ±1). Since in this example —1 is a defective eigenvalue, 
the approximations returned by polyeig have comparable absolute errors 
of order 10^^ which are in accordance with the unstructured backward 
error estimates given in [l^ . 
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2. The accuracy of the remaining approximations is unaffected from the oc- 
currence of near-to-critical eigenvalues and is in accordance with the re- 
sults returned by polyeig. For most non-exceptional eigenvalues, the ac- 
curacy of approximations computed by our method is slightly better. 

3. This kind of behavior is confirmed by many other experiments. Our 
method performs similarly to the QZ for non-exceptional eigenvalues and 
for defective exceptional eigenvalues, but generally worse than QZ and the 
structure-preserving methods [25j for exceptional eigenvalues. 

8. Conclusions and future work 

In this paper we have shown that the Ehrlich-Abcrth method can be used 
for solving palindromic and T-palindromic generalized eigenproblcms. The basic 
idea can be applied to a generic matrix polynomial of any kind; moreover, as 
we have shown in this paper, it is possible to adapt it in order to exploit certain 
structures as the palindromic structure that we have considered here. The 
resulting algorithm is numerically robust and achieves computational efficiency 
by exploiting the rank-structure of the associated linearization in the Dickson 
basis. The algorithm is quite interesting for its potential for parallelization 
on distributed architectures and, moreover, can be easily incorporated in the 
MPSolve package to develop a multiprecision root finder for matrix polynomial 
eigenproblcms. 

There are, however, some issues that still stand in the way of a fully satis- 
factory implementation of our method and are currently under investigation. 

1. The development of an automatic procedure for the selection of starting 
points is important to attain a low operation count due to the prevalence 
and ease of deflation. We have shown that a smart choice could be based 
on a few parameters to be determined from some rough information on 
the spectrum localization. 

2. The proposed algorithm is still inefficient with respect to the size of the 
polynomial coefficients. The preliminary reduction of the linearized prob- 
lem into a Hessenberg-triangular form is the mean to devise a unified 
efficient algorithm for both small and large coefficients. A fast reduction 
algorithm would be incorporated in our implementation. The algorithm 
should be able to exploit the rank structure of the linearization (for large 
degrees), and, at the same time, the inner structure of the quasiseparable 
generators (for large coefficients). 

3. Regarding the accuracy of the method there are still some difficulties in 
the numerical treatment of the critical cases. Our current research is 
focusing on the issue of a structured refinement of the approximations of 
such eigenvalues. 
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